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The conventional theory of homogeneous and heterogeneous nucleation in a supersaturated vapor 
is tested by Monte Carlo simulations of the lattice gas (Ising) model with nearest-neighbor attractive 
interactions on the simple cubic lattice. The theory considers the nucleation process as a slow (quasi- 
static) cluster (droplet) growth over a free energy barrier AF* , constructed in terms of a balance of 
surface and bulk term of a "critical droplet" of radius R* , implying that the rates of droplet growth 
and shrinking essentially balance each other for droplet radius R = R* . For heterogeneous nucleation 
C**) , at surfaces, the barrier is reduced by a factor depending on the contact angle. Using the definition of 

"physical" clusters based on the Fortuin-Kasteleyn mapping, the time-dependence of the cluster size 
distribution is studied for "quenching experiments" in the kinetic Ising model, and the cluster size 
£* where the cluster growth rate changes sign is estimated. These studies of nucleation kinetics are 
compared to studies where the relation between cluster size and supersaturation is estimated from 
' equilibrium simulations of phase coexistence between droplet and vapor in the canonical ensemble. 

The chemical potential is estimated from a lattice version of the Widom particle insertion method. 
For large droplets it is shown that the "physical clusters" have a volume consistent with the estimates 
1 from the lever rule. "Geometrical clusters" (defined such that each site belonging to the cluster is 

occupied and has at least one occupied neighbor site) yield valid results only for temperatures less 
than 60% of the critical temperature, where the cluster shape is non-spherical. We show how the 
'-4_ chemical potential can be used to numerically estimate AF* also for non-spherical cluster shapes. 
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I. INTRODUCTION 



Since the theory of nucleation phenomena was introduced a long time ago [l|43|, the question under which conditions 
the "conventional theory" of nucleation is accurate has been debated (see e.g. [4T-[24||) and this debate continues until 
today. For the simplest case of homogeneous nucleation (by statistical fluctuations in the bulk) of a one-component 
liquid droplet from the vapor, the basic statement of the theory is that under typical conditions nucleation processes 
are rare events, where a free energy barrier AF* very much larger than the thermal energy fcgT is overcome, and 
hence the nucleation rate is given by an Arrhenius law, 



j=uexp(-AF~/k B T). (1) 



Here j is the number of nuclei, i.e. droplets that have much larger radii R than the critical radius R* associated with 
the free energy barrier AF* of the saddle point in configuration space, that are formed per unit volume and unit time; 
w is a kinetic prefactor. Now AF* is estimated from the standard assumption that the formation free energy of a 
droplet of radius R can be written as a sum of a volume term (oc 47rF! 3 /3), and a surface term (oc 4irR 2 ), i.e. 

^ ' 4ttF 3 

; AF(R) = — Afx(p e - Pv ) + 4irR 2 lve . (2) 

• • 6 

>: 

Since the liquid droplet can freely exchange particles with the surrounding vapor, it is natural to describe its ther- 
modynamic potential choosing the chemical potential p, and temperature T as variables, and expand the difference 
in thermodynamic potentials of liquid and vapor at the coexistence curve, Ap = fj, — p coex , Pv and pe denoting the 
densities of the coexisting vapor (v) and liquid (£) phases. According to the capillarity approximation, the curvature 
dependence of the interfacial tension is neglected, is taken for a macroscopic and flat vapor-liquid interface. 
Then the critical radius R* follows from 



dAF(R) 



OR 

and the associated free energy barrier is 



°> R * = A f 1Vl \ ■ ® 

/\fi{Pt - p v ) 



Air 

AF* om = -{R*f lvl . (4) 



However, since typically AF£ om is less than 100 fcnF the critical droplet is a nanoscale object, and thus the treatment 
Eqs. (HI"® is questionable. Experiments (e.g. [25 27]) were not able to yield clear-cut results on the validity of 
Eqs. (HJ)- , and how to improve this simple approach: critical droplets are rare phenomena, typically one observes 
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only the combined effect of nucleation and growth; also the results are often "contaminated" by heterogeneous 
nucleation events due to ions, dust, etc. |28l - l33 |. and since j varies rapidly with the supersaturation, only a small 
window of parameters is suitable for investigation. Therefore this problem has been ver y a ttractive, in principle, for 
the study via computer simulation. However, despite numerous attempts (e.g. fil-fiol [iH l2l| - [23 . [33|V this approach is 
also hampered by two principal difficulties: 

(i) Computer simulations can often only study a small number of decades in time, 0, HH , which in typical cases 
correspond to small barriers (AF* < lOfc^T) rather than the larger ones which are of more interest in the 
context of experiments. 

(ii) On the atomistic scale, it is a difficult and not generally solved problem to decide which particles belong to 
a droplet and which particles belong to its environment; the vapor- liquid interface is diffuse and fluctuating 

[Mill. 

For these reasons, many of the available simulation studies have addressed nucleation in the simplistic Ising (lattice 
gas) model, [l,[l3, (3111,133, 35 56], first of all since it can be very efficiently simulated, and secondly because one can 
define more precisely what is meant by a "cluster". Associating Ising spins oi = +1 at a lattice site i with a particle, 
Ui = —1 with a hole, originally "clusters" were defined as groups of up- spins such that each up-spin in a cluster has 
at least one up-spin as nearest neighbor belonging to the same cluster [331 ] . However, now it is well understood that 
these "geometrical clusters" in general do not have much physical significance |57H6l| : e.g., it is known that there 
exists a line of percolation transitions, where a geometrical cluster of infinite size ap pears, in the phase diagram [57| . 
This percolation transition is irrelevant for statistical thermodynamics of the model |62h65| . 

Based on the work of Fortuin and Kasteleyn (66[ 67] on a correlated bond-percolation model, it is now understood 
that physically relevant clusters in the Ising model should not simply be defined in terms of spins having the same 
orientation and are connected by nearest neighbor bonds, as is the case in the "geometrical clusters" , but in addition 
one has to require the bonds to be "active" : bonds are "active" with probability p 

p(T) = 1 - exp(-2 J/k B T) , (5) 

J being the Ising model exchange constant. 

Due to Eq. ([5]), the "physical clusters" defined in this way are typically smaller than the geometrical clusters, and 
their percolation point can be shown to coincide with the critical point (59j - [6l| . A geometrical cluster hence can 
contain several physical clusters. Note that to apply Eq. ([5]), random numbers are used, and hence physical clusters 
are not deterministically defined from the spin configuration, but rather have some stochastic character. This presents 
a slight difficulty in using physical clusters in the study of cluster dynamics. 

While Eq. ([SJ has been used in the context of simulations of critical phenomena in the Ising model, applying very 
efficient Swendsen-Wang [6(| and Wolff [g8[ simulation algorithms, this result has almost always been ignored in the 
context of simulations of nucleation phenomena [IH, [22|, |47T - |52| . While it is allright to ignore the difference between 
geometrical and physical clusters in the limit T — > (obviously p(T) — > 1 then, all bonds becoming active), this is 
completely inappropriate at higher temperatures. 

The present work hence reconsiders this problem, studying both dynamical aspects of nucleation in the framework 
of the kinetic Ising model [6^, [7(3] (without conservation laws), and the static properties of large droplets, applying 
the definition of "physical clusters" based on Eq. ([5]) throughout. For comparison, we shall also occasionally use the 
"geometrical" cluster definition, to demonstrate that misleading conclusions would actually result in practice, for the 
temperatures that are commonly studied. The study will be generalized to Ising systems with free surfaces, where a 
boundary field Hi acts [H, . First of all, in this way also a systematic investigation of heterogeneous nucleation at 
planar walls becomes feasible; secondly, due to the reduction of the barrier AF^ ot in comparison to Ai 7 j^ om ; nucleation 
for reasonably large values of R* becomes accessible to study. 

In Sec. |TT1 we consider the equilibrium of the lattice gas model for p v < p < pi in systems inaLxLxL geometry 
with periodic boundary conditions, to show that physical clusters do occupy precisely the volume predicted by the 
lever rule analysis, [U [23|, [24], Hfl as they should when the thermodynamic limit is approached. We present evidence 
that physical clusters are correctly identified by both the lever rule method and the approach based on the "atomistic" 
identification of clusters based on Eq. ([5]) at all temperatures, from zero temperature up to the critical temperature 
T c . In contrast, Eq.Q, which implies a spherical droplet sha pe, is found to work only at temperatures distinctly 
above the interface roughening transition temperature Tr |7ll [72j, even for very large radii R. We attribute these 
discrepancies to the fact that due to the anisotropy of the interface tension for our lattice model pronounced deviations 
of the average droplet shape from a sphere occur [73l - [76| . presenting data on the shape of large droplets. In Sec. IIII1 
we describe our results on the dynamics of the droplet size distribution and on the attempt to find R* from the size 
where growth and shrinking processes of clusters are balanced. This study is also carried out for systems with a free 
surface, for which the contact angles for various values of the surface field have been estimated previously [55L [56| . 
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since in this case much lower barriers (for large clusters) result, which is crucial for making this study feasible with 
manageable effort. However, the radii R* predicted from this analysis of kinetics show slight deviations from the radii 
R* predicted from AF^ ot . Possible reasons for this discrepancy will be discussed. Finally, Sec. ITVl summarizes our 
conclusions. 



II. MICROSCOPICALLY DEFINED CLUSTERS VERSUS MACROSCOPIC DOMAINS IN THERMAL 

EQUILIBRIUM 

As is obvious from Eq. ([2]) and the reasoning behind it, this approach is adequate when one deals with the description 
of macroscopically large domains in equilibrium with a surrounding bulk phase. However, one needs to find an 
extension of the concept that can be applied also to nanoscopically small droplets, "clusters" in the lattice gas model 
that contain perhaps only of the order of 100 fluid particles. In this section, we want to confirm the idea that one must 
use the concept of "physical clusters" based on Eq. (0 for this purpose, rather than the "geometrical clusters" that are 
so widely used when the lattice gas model is used to test nucleation theory concepts. While the geometrical clusters 
are appropriate if one works at extremely low temperatures where the clusters basically have the shape of small cubes 
[52j . this region clearly is inappropriate when one has the application for vapor-to-liquid nucleation in mind, where 
droplets are spherical, and their interfaces are rough and fluctuating rather than smooth planar facets. In fact, many 
studies of nucleation in the lattice gas model have been made in d = 3 dimensions at temperatures near T/T c = 0.6 
or thereabout; given the fact that the interfacial roughening transition of the Ising model on the simple cubic lattice 
is known to occur at about [77], III] T R /T C « 0.544, i.e. (note k B T c /J = 4.51154 79]) k B T R /J w 2.44, it is clear that 
temperatures much closer to T c must be studied to render the assumption of a spherical droplet shape accurate. In 
fact, this assumption of a spherical droplet shape is accurate when the difference between the interfacial stiffness [8(| 



and the interfacial free energy becomes negligibly small. Numerical studies of Hasenbusch and Pinn [8l| indicate that 
this is only the case for ksT / J > 3.9. As a consequence, it is clear that most of the existing studies of nucleation 
phenomena in the Ising model, that were based on geometrically defined clusters, and had to be done at much lower 
temperatures, are inconclusive: the deviation of the average droplet shape from a sphere enhances the surface term 
in Eq. ([2]) ; but the fluctuation corrections discovered for small droplets by the "lever rule method" [U [23|, [24], HH 
show that y v e(R) < j v i(oo) for small R and hence the surface term in Eq. ([2} is decreased. Thus, it hardly can be a 
surprise that some of the studies concluded that the nucleation barriers predicted by classical nucleation theory and 
the capillarity approximation (that ignores the i?-dependence of y v g(R)) are too high, and others concluded they are 
too low, or even reported good agreement. We take the latter finding as indication that the two opposing effects have 
accidentally more or less canceled each other. 

This problem is the motivation for the present section, which attempts to show that "physical clusters" based on 
Eq. ([5]) are appropriate to identify clusters in the lattice gas model, irrespective of temperature and cluster size, and 
are equivalent to the droplets of the "lever rule method" , for large enough droplets. 

Fig. [1] recalls this approach: one samples for a system of volume V = L x L x L the effective thermodynamic 
potential /j, (T, p) per lattice site as a function of density p, 

f L (T, p) = [F(N, V, T) - F(N = Vpi, V, T))/V , (6) 

where N is the number of occupied lattice sites (pi = (1 + <7j)/2 = +1 with Oi = ±1 the spin variable at lattice 
site i). Since phase coexistence between bulk liquid (at density pi) and vapor (at density p v ) occurs at a chemical 
potential p C oex that corresponds to the "field" (in magnetic notation) H = 0, pi and p v are simply related to the 
spontaneous magnetization of the Ising ferromagnet as p£ = (1 + m sp )/2, p v = (1 — m sp )/2, and H translates into p 
via H = (p, — /! CO cx)/2 = A/i/2. The accurate sampling of fL,(T, p) for large L is a nontrivial task, it requires the use 
of advanced methods such as "multicanonical Monte Carlo" [93[ or "Wang Landau sampling" [9J, [95| or "successive 
umbrella sampling" [9(3, HH, see [H, [HJ for background on such techniques. From Eq. ©, one defines a chemical 
potential function as a derivative, 



p{N,V,T) 



dF(N, V, T) 



(7) 

V,T 



dN 

Mcocx ■ (8) 

One recognizes that the isotherms A/iz,(T, p) vs. p exhibit a loop: the homogeneous vapor remains stable also for 
some region where p > /i CO cx, until a peak occurs, which indicates the "droplet evaporation/condensation transition" 
[H], [111 HIJ : in the first regime where A/i^(T, p) decreases with p, a (more or less spherical or cubical) droplet coexists 
with surrounding vapor. Here, we are not interested in the further transitions that one can re cog nize from this curve, 
where the droplet changes shape from spherical to cylindrical, or to a slab configuration, etc. j2ll.l23l.l24l]. Instead, we 
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FIG. 1: Illustration of the numerical procedure for determining the density triplets p' , p and p" and the associated values 
/z,(T, p'), }l(T,p) and /l(T, p") for a given choice of T and L, according to the "lever rule method". The density p must be 
chosen such that it is not too close to the peak of the curve A/ii/fcsT versus p, which is due to the "droplet evaporation- 
condensation transition", because all microstates that are sampled over should contain a droplet in the system. It must also 
be chosen not too close to the first kink in the curve Ap^/feaT versus p, which is due to the transition of the droplet to a 
cylindrical shape (stabilized by the periodic boundary condition in the direction of the cylinders). Apart from these constraints, 
the choice of p is arbitrary, and studying different choices of p provides useful consistency checks on the results. The densities p' 
and p" then can be read off from the curve Ap^/ksT versus p, since it is required that Ap£,(T,p') = Apl(T,p) = Ap^{T, p"). 
The corresponding values of the free energy densities then can be read off from the lower part of the figure, and using Eqs. ([9]), 
pop with iV cxc = 0, the surface free energy F$ is extracted. The actual data shown here refer to the case IcbT/J — 3.0 and 
L — 20. 



emphasize the key idea of the "lever rule method" : one can identify a range of choices for the chemical potential p 
where three states of the finite system can exist in equilibrium with the same chemical potential, namely a homogeneous 
vapor at density p' > p v , a homogeneous liquid at density p" > pi, and a state where two-phase coexistence between 
the droplet and surrounding vapor occurs. Since the vapor in this case exists at the same chemical potential as the 
pure vapor, it must be of the same physical nature as the state with density p' , and similarly, the liquid in the droplet 
can be identified with the liquid at p" . Making now use of the fact that for large enough systems a system can be 
suitably decomposed into independent subsystems, we write for the free energy, with V" the volume taken by the 
droplet, V = V - V", 

Vf L (T, p) = V'f L (T, p') + V"f L (T, p") + ^ s , (9) 

where the free energy densities /l(T, p') and /t(T, p") are explicitly known, and also /l(T, p) is known: thus, when 
the droplet volume V" is known, the surface free energy F$ of the droplet, which is defined via Eq. ((9]), is determined. 
A similar decomposition can readily be written down for the particle number, 

N = Vp = V'p 1 + V"p" + N exc , (10) 

where we have allowed for an excess number N exc of particles, to be associated with the interface. If we consider the 
definition of an "equimolar dividing surface" 34], N CKC — 0, and then reading off p' and p" from the construction in 
Fig. [T] we see that Eq. (jTT)]) readily yields V and V" for the considered density p, and via Eq. © we can immediately 
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extract Fs from the data. Note that these arguments do not invoke the assumption that the dividing surface needs 
to be a sphere. If one makes the assumption, V" — 4irR 3 /3, and then one can write also Fg — 4ttR 2 ^ v ((R). Thus, 
it is assumed that all interactions of particles inside the droplet (volume region V") with particles inside the vapor 
(volume region V') are restricted to the interfacial region, and hence can be accounted for by their contribution to 
the surface free energy F$- 

(a) (b) 
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FIG. 2: (Color online) (a) Plot of p'{T, Ap) versus Ap/ksT, for several temperatures ksT/J as indicated. The symbols 
represent data recorded in the grand-canonical ensemble of the lattice gas, while the curves are obtained in the canonical 
ensemble, recording Ap — Ap(T, p) via the lattice version of the Widom particle insertion method [5oT. l87l. l8Sj . The chosen 
lattice size was L = 60. (b) Same as (a), but for p" (T, Ap). Both methods agree perfectly. 



However, the method defined via Eqs. ©, (fTUj) . and illustrated in Fig. Q] becomes difficult to apply for very large 
droplets, because the sampling of /z,(T, p) then becomes unreliable or would require an unaffordable effort. The 
method is also difficult to apply for small droplets, because then one must use relatively small simulation boxes to 
ensure the stability of the inhomogeneous state where a droplet coexists with surrounding vapor [37| . Thus, it is 
very desirable to complement the approach by a more "microscopic" identification of droplets, and this is possible 
via Eq. ([5]). In particular, it has been shown that apart from finite size effects (see e.g. [61() that for L — > oo 
the spontaneous magnetization of the Ising ferromagnet m sp coincides with the percolation probability P, which 
is defined [86( as the fraction of sites belonging to the largest "physical cluster" in the system. When we hence 
analyze a configuration at a density p, where (cf. Fig. Q]) a large cluster is present in the system, using Eq. (|5|) to 
define clusters the largest cluster will include I sites, which we hence can associate with its (total) magnetization 
M = m"V" where m" is then the magnetization per site (note that the Ising magnet/lattice gas isomorphism implies 
that p" = (I + to")/2). As a consequence, we can obtain the droplet volume V" from a "measurement" of the average 
size (£) of the largest cluster in the system via 

v = 4 ■ (ii) 

to" 

Note that this volume in general differs from the volume of a geometrical cluster: If (l gcom ) counts all occupied sites 
belonging to the geometrical cluster, and noting that the density in a large geometrical cluster is just the bulk density, 
namely p" = (1 + to")/2, the volume taken by the geometrical cluster is given by 

yll _ (Igcom) _ 2(l gcom ) 

"ge°m- pll ~ 1 + m » ' 1 > 

In the limit L — > oo, where also (£} and V" get macroscopically large, m" tends to m sp , while for finite L it is clear 
that to" slightly exceeds m sp (and p" exceeds pi, see Fig. Q]). However, recording the relation p = p(T,Ap,) [or the 
equivalent relation to = to(T, H) of the Ising ferromagnet] very precisely is an easy task, see Fig. [21 since both states 
at densities p',p" (Fig. [T]) are homogeneous, not affected by heterophase fluctuations, and since the temperatures 
studied are still well below T c , statistical fluctuations are small, and finite size effects arc negligible. Fig. [2] presents 
representative results for both p'{T Ap) and p"(T, A/x) versus Ap,. Note we also have used the lattice version of the 
Widom particle insertion method [HH [87l . f88j j to record the inverse function A/i = A/i(T, p) from simulations in the 
canonical ensemble, where p was chosen as the independent control variable. The perfect agreement between both 
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approaches not only serves as a test of the accuracy and correctness of our numerical procedures, but also shows 
that for the chosen temperatures and linear dimensions finite size effects on states in "pure" phases are completely 
negligible, since finite size effects are known to differ in the two ensembles (H!, [H[, but are not detected here at all. 
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FIG. 3: (Color online) Plot of the ratio V"/V/ r ' of the droplet volume as obtained from the "Coniglio-Klein-Swendsen-Wang" 
[59l. |60| cluster definition, using Eq. (J5]) to define physical clusters, and their volume V" (Eq. HI])), and from the lever rule 
(Eq. (|13|) h as a function of 1/L, for four temperatures: fcsT/J = 2.0 (a), 2.6 (b), 3.0 (c) and 4.0 (d). Various densities are 
included for each temperature, as indicated. 



To record V" as defined in Eq. (|11[) from the simulations, we performed simulations in the canonical ensemble, 
choosing various values of p, and equilibrate a large droplet coexisting with surrounding vapor. The initial state 
is then chosen putting a droplet with the size predicted by the lever rule (and density p = (1 + m sp )/2) into 
the simulation box, which then is carefully equilibrated. The standard method to simulate the Ising model with 
conserved magnetization (which corresponds to the lattice gas model in the canonical NVT ensemble) is the "spin 
exchange algorithm" [83, 84]. However, the standard nearest neighbor exchange method implies that any local excess 
of magnetization (or density, respectively) can only relax diffusively, and the resulting "hydrodynamic slowing down" 
[83l [84| hampers the fast approach towards thermal equilibrium. We thus used instead a single spin-flip algorithm in 
which the total magnetization is restricted to two neighboring values {M — 2, M}: So if the system has magnetization 
M, flipping a down-spin (which would mean a transition M — > M + 2) is automatically rejected, and if the system is 
in the state M — 2, flipping of up-spins is forbidden. For large systems (L — > 00) any corrections to a strictly canonical 
simulation at a magnetization per spin m = M/L 3 are of order 1/L 3 and hence negligible. 

Choosing very large systems (up to L = 160, rather than L = 20 as used in Fig. [IJ the ratio of V" as found from 
Eq. ((TT|) and VJ" from the lever rule (with the assumption -/V exc = 0), i.e. 



v" = v p - p ' 

lr P"~P" 



(13) 
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is plotted vs. 1/L in Fig. [3] for several temperatures and various densities p. These data show that V" /V{' — > 1 as 
L — > oo, irrespective of temperature and density (in the density region where a droplet not affected by the periodic 
boundary conditions is present, as explained in Fig. H}. The fact that V" /V/l extrapolates to unity for L — » oo not 
precisely, but only within some error, is due to the fact that statistical errors affect both the estimation of {£) and of 
V(' r (via errors in the estimation of A/i and hence p'). As expected, using the proper definition of physical clusters, 
one can work at arbitrary temperatures, both at T < Tr (a), T slightly above Tr (b) or T rather close to T c (d). 
While in cases (a) the simple geometrical cluster definition would also work, since essentially all bonds are "active" 
(p(T) is almost unity {Eq. ©}), and the non-spherical shape of the clusters does not matter in this context. But the 
geometrical cluster definition would clearly break down in case (d) due to the proximity of the percolation transition 
that occurs for geometrical clusters at a density not much larger than those included in Fig. [3Jd) [13, Ifill, [65J . On 
the other hand, we note from the fact that there always occurs asymptotically in the relation V" /VI" a correction of 
order 1/L, that for physical clusters the assumption N cxc — 0, that is often (but not always [24]) made in the lever 
rule method, does not hold: i.e., when we assume that N cxc is proportional to the surface area of the droplet, we can 
write 



N cxc = C(V") 2/3 Pc 



(14) 



where C is a geometrical factor (C = (367T) 1 / 3 for a spherical droplet, C — 6 for a cube), and p cxc is an excess density 
of the particles due to the interface of the droplet. Using V/l from Eq. ([TBI as a first-order estimate in Eq. (Ti"4]) . one 
readily finds that the term N cxc /V in Eq. (fTCTj) yields a 1/L correction, 
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P ' 



V L \p" 
and hence Eq. (fTU)) would yield instead of Eq. (pi 



P' 



2/3 



(15) 



V" r (corrected) 



V," 



c 

= i — 

L 
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2/3 



= 1 - 



P - P J P- P 
Cpexc^/i _ /\-2/3( 



L 



(16) 



-(p"-pT 2/3 (p-pr 1/3 



which is qualitatively in accord with Fig. G3 In order to test Eq. (fT6| and obtain estimates for the temperature 
dependence of p oxc , Fig.H]plots our numerical results for the ratios V" /V t " versus (p" — p')~ 2 ^ 3 (p — p')^ 1 ^ /L. We see 
a very good data collapse at straight lines going through unity at the ordinate within numerical error at all studied 
temperatures. 
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FIG. 4: (Color online) Plot of V" /Vi" versus (p" — p')~ ' (p — p')~ ' /L at low temperatures (a) and at temperatures closer to 
T c (b), including all densities p that were analyzed. Straight line fits to the data are included. Note that the different symbols 
correspond to different choices of p in Fig. [3] 



Since the data in Fig. [3] suggest that for "physical droplets" in the lattice gas model an appreciable "interfacial 
adsorption" (as expressed in p oxc in Eq. (I14[) ) occurs, it is of interest to not only study the average volume of the 



8 



(a) 



(b) 




(d) 




35 
30 
25 
20 
15 
10 

5 - 





30 40 
Radius R 



T=3 
T=4 
T=4.3 



-l I + + + + + I I + + 



+++*+ ■ 



+ + l + 



Radius I 



2.8 3 3.2 3.4 3.6 

logarithm of droplet radius log(R) 



3.8 



FIG. 5: (Color online) Radial density profiles p(R) plotted vs. 7? at k B T/J = 3.0 (a), T = 4.0 (b) and k B T / J = 4.3 (c), for a 
box of size Lx Lx L with L = 60, 80, 100, 130, 160 (from left to right). The vapor density and the liquid density are independent 
of L, as expected. Vertical straight lines indicate the radii that follow (assuming R = (3V" /4ir) 1 ' 3 , i.e. a spherical shape) from 



Eq. ()lip with V" — (l)/m" or from the simple geometrical cluster definition, Eq. (|12|) . where Vg eom = 2(£ gcom )/(l + m") or 
from the lever rule Eq. ()13|) . A horizontal line at p — 0.5 is also drawn, since for R — > oo the density profile should become 
symmetric with respect to this line, due to the spin reversal symmetry of the Ising model, and then the correct cluster definition 
must yield a cluster radius compatible with this inflection point of the profile. All data are obtained from averages over several 
hundred independent droplet configurations. The dotted vertical lines are the radii predicted from the lever rule. Part (d) gives 
a plot of the average squared interfacial width (w 2 ) vs. InR (in this plot, R is taken from the condition p(R) = 0.5). 



droplets (£) but also their radial density profile (Fig.[S]). It is seen that for radii R which are in the range from 10 to 
20 lattice constants (corresponding to droplet volumes in the range from about 4500 to 36000, so these are already 
clusters of a mesoscopic size, with a huge nucleation barrier far beyond observation in a simulation of nucleation 
events or in experiments) the profiles are very broad, and their width increases slightly with increasing droplet size. 
For comparison, the prediction for the cluster radius resulting from the widely used standard geometrical definition of 
clusters is also included: it happens that this geometrical radius is still fairly close to the correct radius at ksT/J = 3.0. 
Of course, closer to T c the geometrical cluster definition yields completely unreasonable results, due to the onset of 
percolation phenomena, and for fc^T/ J = 4.0 and 4.3, the geometrical radii are indeed unreasonably large. 

It is interesting to note (Fig. [SJi) that the width of the interfacial profile increases with increasing R. This phe- 
nomenon is well-known for planar interfaces and attributed to capillary waves. Since for a large droplet the surface 
is locally planar, most of the capillary wave spectrum is not affected by the interface curvature. Thus we may, as a 
first approximation, take over the result for the broadening of a planar interface of linear dimension L, replacing L 
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by the droplet radius R [98l-fl00|] 

where wo is the "intrinsic width" of the interfacial profile, 7 the "interfacial stiffness" (note that a factor 
is absorbed in its definition) and A m i n a short wavelength cutoff, whose precise value is not known. Near T c the 
interfacial stiffness coincides with the interfacial free energy, while 7 — ¥ 00 at T — > Tr, and then the capillary wave 
broadening disappears. If one accepts the above formula, and uses the data of Hasenbusch and Pinn [8lj ] . one predicts 
for the slope [47] _1 sa 2.5 (k B T/J = 4.0) or 7.98 (k B T/J = 4.3), respectively. The actually observed slopes of the 
lni? term are actually somewhat smaller, namely about 2.13 at fc^T/J = 4.0 and 5.83 at /cgT/J = 4.3, but of the 
same order of magnitude. 

We deliberately do not discuss the radial droplet density profiles for ksT/J = 2.0 and ksT/J = 2.6, however, since 
at these temperatures the droplet shape shows distinct deviations from the spherical shape. This can be checked 
directly by recording contours of constant density in slices of width = 1 (taking into account 3 lattice planes through 
the droplet's center of mass, parallel to the xy plane, the xz plane and the yz plane, respectively). Averaging these 
density contours over several hundred statistically independent observations the plots shown in Fig. [5] are obtained. 
They are all taken at the same fixed density p = 0.33. Note that for this density, the stable state would be a cylindrical 
droplet (stabilized by the periodic boundary conditions) but for such large systems (L = 1 00) the compact droplet 
shapes shown here (chosen by an appropriate initial condition) are always perfectly metastable. At fc^T '/ J — 2.0 
the cubic symmetry of these cross sections through the droplet is evident (although the presence of facets parallel 
to the planes x = or y = is not evident, due to finite-size effects at the points where the facets join the round 
sections, which replace the sharp edges of the cubes at nonzero temperature). At ksT/J = 2.6, where we exceed the 
roughening temperature slightly, there clearly occur no longer any facets, but the density contours in Fig. [5] are still 
distinctly non-circular: the diameter in diagonal direction clearly is about 7% larger than in the lattice directions. 
Even at ksT/ J = 3.0, we find an enhancement of the diameter in diagonal direction of about 3%. At fc^T/J = 4.0 to 
4.3, however, no longer any statistically significant deviation from spherical droplet shapes (and hence circular shape 
of the density contours in the cross sections) can be detected. Note that for ksT/J = 4.3 (case (f)) the geometrical 
cluster definition would not be applicable due to the proximity of the percolation transition of geometrical clusters. 
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FIG. 6: (Color online) Contours of constant density p(x,y) — pi in the planes x — 0, y — and z — cutting through 
droplets situated at the origin in a box of dimension L = 100 with periodic boundary conditions in every direction, for the 
cases k B T/J = 1.0 (a); k B T / J = 2.0 (b); k B T/J = 2.6 (c); k B T / J = 3.0 (d); k B T / J = 4.0 (e); k B T / J = 4.3 (f). The 
density is always fixed at p = 0.33. In addition to the color-coded average density per site, the contours for three densities 
Pi = (pviT) + 0.5)/2, 0.5, (pe{T) + 0.5)/2 are shown. 
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FIG. 7: (Color online) Droplet volume V" relative to the lever rule estimate V{ r plotted versus 1/L for four temperatures 
UbT/J = 2.6 (a), 3.0 (b), 4.0 (c) and 4.3 (d). Three definitions of the droplet volume are compared to each other: the 
"geometrical cluster" definition Vg" om = 2{£ gC om)/(l + m"), the "physical cluster" definition V" = {£)/m" {Eq. (|lll) } and 
the result of classical nucleation theory V" = 47r(iZ*) 3 /3 where R* is computed from Q, using (A/i) as "measured" in the 
simulation by the lattice version of the Widom particle insertion method. At each temperature, several densities were used: 
p = 0.05, 0.06, . . . , 0.12 (a), 0.07, 0.08, . . . , 0.12 (b), 0.17, 0.18, . . . , 0.20 (c), 0.26, 0.27, 0.28 (d). 



As a final part of our analysis of static properties of physical droplets in the Ising model, we exploit the fact that the 
chemical potential p, can be "measured" by the lattice version of the Widom particle insertion method [5^, H3, HI| also 
in the state when the system is inhomogeneous, e.g. for the case of interest when a droplet coexists with surrounding 
vapor. Actually, already in [55j it was shown that p, actually stays spatially constant in such a situation. However, 
while the estimation of p from Eq. (|7J| requires a very accurate estimation of the free energy density T) {Eq. ([5]), 
Fig- D]}, an( i m practice this works only for not so large linear dimension L (such as L = 20 in Fig. [IJ, the estimation 
of p, from the particle insertion method still works for volumes that are orders of magnitude larger. As a consequence, 
we can use Eq. ([3]) to test whether or not the actual droplet volumes V" are compatible with conventional nucleation 
theory for large droplets (assuming that the droplets are spherical) so R* — (3V" /4-it) 1 / 3 holds for a critical droplet. 

Fig. [7] presents a plot of the cluster volume V" versus inverse linear dimension for several densities, at the temper- 
atures ksT/J = 2.6, 3.0, 4.0 and 4.3, using the classical nucleation theory prediction Vq NT = 47r(i?*)'V3 with R* 
given by Eq. (J3j) for comparison. For this purposes j V £ is taken from the work of Hasenbusch and Pinn (8ll . 82], and 
so there occur no unknown parameters whatsoever. The geometrical cluster volume V^" om {Eq. (|12p | is close to the 
estimate based on the Coniglio-Klein-Swendsen-Wang "physical cluster" -definition, Eq. (fTTj) . at fc^T/J = 2.6 and 
fc^T/J = 3.0, while at ksT/J = 4.0 (and higher) the deviations become appreciable: V^" om then is systematically too 
high (in comparison with all other estimates, including VJ" {Eq. (fl"3)) }. which is used as a convenient normalization). 
It is interesting to observe that the classical nucleation theory estimates based on Eq. © are systematically too low 
for T = 2.6 and T — 3.0, while for T > 4.0 Eqs. (fT3")) are found to be in very good agreement. This discrepancy 
at the relatively low temperatures comes from the fact that using Eq. ([3]) for the estimation of V" we imply that 
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the cluster volume is spherical and hence we underestimate the surface area: the geometrical factor C (T) introduced 
in Eq. ([T4")) increases from about 4.836 for the spherical shape near T c up to 6.0 for the cube, as the temperature is 
lowered. Since the non-spherical droplet shapes at k B T / J = 2.6 and 3.0 have regions of rather large curvature near 
the parts of the droplet where at low temperature the edges of the cube will appear, the estimation of p then yields a 
too low radius R* . The observation that the anisotropy of surface tension in the lattice gas model becomes noticeable 
for k B T/J < 4.0 is consistent with the findings of Hasenbusch and Pinn [8ll j . 

If all the methods to define clusters were correct, in the thermodynamic limit all data should extrapolate to 
V" /V(l — 1 for L — > oo, since in this limit the lever rule {Eq. (1131) } is trivially true with p' — p v and p" — p£. The 
method based on the definition of physical clusters, Eq. (fTTj) . indeed is nicely compatible with this expectation at 
all temperatures; although it is somewhat unsatisfactory (and unexpected) that at finite L there occurs a surface 
correction due to the surface excess 7V exc noted in Eq. (TPfj) (and discussed above). However, it is clear that the two 
other methods do not give results that are correct for L — > oo in general: while the method based on the geometric 
cluster definition still gives essentially correct results at k B T/J = 2.0 (not shown here), where the distinction between 
"geometrical" and "physical" clusters is irrelevant, for temperatures T > Tr, the volume of geometrical clusters is 
systematically too large. At k B T/J = 4.3, the error is as large as 60 to 80% even asymptotically, and for the cluster 
sizes that were actually studied the overestimation actually is by a factor two to five (Fig.[7JI)! In view of the fact that 
the geometric cluster definition must break down due to the percolation transition [57] , this failure is not unexpected, 
but we are not aware that it ever has been quantified previously. In the regime from fc^T/ J = 2.6 to 4.0, the error of 
the geometric cluster volume raises from a few percent to 15 to 40%. 

The results obtained from the classical nucleation theory via the "measurement" of the supersaturation Ap, on the 
other hand, yield essentially the correct result at high temperatures (k B T/J > 4.0), for all linear dimensions studied, 
since the data are essentially independent of L, and hence R*, in the considered range. But it is remarkable that 
for k B T/J = 3.0 we find Vg NT /V£ ~ 0.94 (Fig. [7b) and for k B T/J = 2.6 we find V^ T /V{1 w 0.89 (Fig. [7k). This 
discrepancy becomes worse at lower temperatures (e.g. Vq NT /V^" « 0.68 at k B T / J = 2.0 [not shown]), and obviously 
this discrepancy must be attributed to the orientation dependence of the intcrfacial free energy and the resulting 
non-spherical droplet shapes (Fig. [5]). 

At first sight, the result that Vq^/V^ is independent of R* seems to be at variance with the finding of a curvature- 
dependent surface tension ^ v i{R) due to Winter et al. [2l|, [H, [24l [H, [H, [96| . In fact, evidence was provided that 

lAR) TTWiW ( ] 

where / is a length proportional to the correlation length in the bulk 96j . Note that due to the spin reversal symmetry 
of the Ising model one can show (9l| that a Tolman correction (oc 1/R) must be absent for R — > oo. However, from 
Eqs. @ and (fT5)l . it is straightforward to show that 

R Ap{p e - p v ) = 27(00' 



(l + 2(Z/i?*) 2 ) 2 (19) 
= 2 7 (oo) [l-4(l/R*) 4 + 0((l/R*f)] . 

Hence for R* J> I it is clear that the result for R* is still given by Eq. ([3]), and for k B T/J < 4.0 we are safely in 
this regime, since the correlation length then does not yet exceed the lattice spacing. So the deviations of the ratio 
^c'nt/^" that are seen in Fig. [7J can be attributed fully to the deviation of the droplet shape from a perfect sphere, 
caused by the anisotropy of the interfacial free energy. In Fig. [8l we now present the ratio of the intercepts V£ /Vq NT 
for R* — > 00 as a function of temperature, since we know that {Eq. ©} V" C * NT = (47r/3)7^(oo)/[A/x(# - p v )} 3 , 
while Eq. dU) yielded AF* om /k B T = (367r) 1 / 3 7 t ,^(oo)/3(Vg NT ) 2 / 3 = m coex HV£ NT for a spherical droplet. At T = 0, 
however, the drople t is a perfect cube, and for intermediate temperatures, its shape (for V* 00) is given by the 
Wulff construction [l0l| (and hence not explicitly known). However, for large droplet volume V we can write in 
general 

AF(V) = -2m cocx HV + ~/AV- 2/3 V 2/3 , V -> 00 (20) 

when we have assumed that the droplets of different linear dimension R (for R — > 00) have the same shape at fixed 
temperature, so we can write V = VR? for the droplet volume (V is then formally the volume for R = 1) and the 
surface area is A = AR 2 . For instance, for a sphere we have V^phcrc = 47r/3 and A sp h ere = 4tt, and for the cube 
Kubc = 1 and A — 6. Minimizing AF(V) with respect to V yields 

(y*)V3 = I_Z AV- 2 ' 3 (21) 

3 Tn C n C vH 
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for a general shape, which is in between sphere and cube. The barrier then can be written as 



AF* 



1 



7 



27 (m coex H) 



-A A V~ 



,HV* 



(22) 



Using now the fact that by choosing a particular large droplet volume in our simulation, H is automatically fixed for 
any volume V* due to the thermal equilibrium situation constructed in our simulation. So it makes sense to estimate 
the ratio of barriers as 



AF* 

afT 



v,: 



CNT 



(23) 



CNT 



For T = 0, we know that 7 is again the interface tension of the planar surface, also used in the classical nucleation theory 
for the spherical surface. Hence when we write the ratio of AF*/AF^ NT , using Eq. (f22|) . the term j 3 /(27(m cocx H) 2 ) 
cancels, 
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should be reached for T — > 0, while 
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(Fig. E]) are compatible with this expectation. 



(24) 



1 as T — » T c . Our numerical results 
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FIG. 8: (Color online) Temperature variation of Vj"/Vcnt- The locations of the roughening temperature Tr and the critical 
temperature T c are indicated by dotted vertical lines. As stated in the main text, the values of 7„( are taken from [8^ |. 
using Monte Carlo results for temperatures above fcsT/J = 2.0 and a low temperature series expansion up to 17th order for 
temperatures below ksT j J = 2.0. Note that for temperatures smaller than ksT/J — 1.4, the magnetization at coexistence 
TOcogx is almost indistinguishable from its saturation value m CO cx = 1, and then our implementation of the Widom particle 
insertion method cannot be applied to the chosen system sizes. 
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III. TIME EVOLUTION OF THE DROPLET SIZE DISTRIBUTION AND DROPLET GROWTH RATES 



We now consider time-dependence of the system where we start the system at time t — in an equilibrated state 
in the vapor at /i = fi coox , but switch on at time t = a chemical potential \x > fJ, cocx (or equivalently, a positive 
magnetic field H in the notation of Ising ferromagnet) at which the liquid is the stable phase. As mentioned in the 
introduction, we consider heterogeneous in addition to homogeneous nucleation, choosing a L x L x D geometry with 
two walls at z = 1 and z = D, choosing surface fields Hi, Hp such that Hp favors the vapor but Hi (acting at the 
wall at z = 1) favors the liquid. The reason for this choice is, that by proper choice of Hi one can adjust the contact 
angle 8 at which sessile wall-attached macroscopic droplets can occur. The barrier against heterogeneous nucleation 
Ai 7 ^ is predicted to be very much reduced, in comparison to the barrier AF£ om against homogeneous nucleation, if 



the contact angle is small, since [H, [1^] 



AF h * ct = AF h * om /(0) 



f{6) = (1 - cos 9) 2 (2 + cos 0) /4 . 



(25) 



As shown with the "lever rule" method [55J, [56| indeed rather large wall attached droplets can be simulated in 
equilibrium with supersaturated vapor which have barriers of order Ai^ ot « lOfcsT or so only, for suitable choices of 
Hi and H, and so a comparison with kinetic studies then seems reachable, and varying Hi over some range provides 
an additional variable to test the theory. 
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FIG. 9: (Color online) Mean lifetime of metastable states inlxLxfl Ising systems at ksT/J = 3.0 (a) and fesT '/ J = 4.0 (b) 
as a function of the chemical potential \ijksT, for various choices of Hi as indicated. This lifetime was measured as the first 
time, when the time-dependent magnetization m(t) becomes positive, as described in more detail in the main text. L — 60 and 
D = 30, with Hd/J = —0.9 throughout, and periodic boundary conditions in x and y directions only. The average lifetime 
decreases with higher surface fields Hi and highter chemical potentials A^i/fcgT. 



As a first step, preliminary runs were performed with the single spin flip Metropolis algorithm [83l |84| monitoring 
the average lifetime of the metastable vapor. This was done using a large sample (10 5 ) of equilibrated initial states 
at H = 0, where the considered field H (or chemical potential fi — ^ COG x, respectively) was then switched on and the 
time recorded when the (initially negative) magnetization reaches the value m = for the first time. Fig. |H] shows 
estimates for the resulting mean first passage times for a range of choices of Hi as a function of the field. When this 
"lifetime" of the state with m < (i.e., vapor) does not exceed 10 4 , the system is rather unstable, nucleation occurs 
fast and is followed by fast domain growth as well; such fast decays of unstable systems are not suitable for tests of 
nucleation theory. It is seen, that in a rather narrow interval of fields H (for each value of Hi) the lifetime increases 
from 10 to 10 6 . Such parameter combinations (H,Hi) will be studied in the following only; if we would study cases 
where the lifetime is significantly larger than 10 6 , no critical droplet would be formed during affordable simulation 
times. 

Fig. 1101 then shows typical time evolutions of the size distribution of ng(t) of physical clusters of size £, normalizing 
them by the equilibrium cluster concentration n° e q for H = 0. Here n c e q is defined as the average number of physical 
clusters per lattice site, and we have the sum rule 

N 

p = E^ q ( 26 ) 
t=i 
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FIG. 10: (Color online) Time dependence of the cluster concentration ratio ni(i)/n i l q in equilibrium at H = 0, for the case 
k B T/J = 3, Hi/ J = 0.4, and H/J = 0.15 (a) and 0.17 (b). Curves show the choices I = 40,50, ... , 120 (from bottom to top). 



since every site occupied by a particle must be part of some cluster. Of course, here we are only interested in not 
too small clusters, and hence Fig. [TU] focuses on clusters with £ > 40. In the time evolution of ni(t)/n e e q we recognize 
three regimes: for times of order t < 100, ni{t)/n cq is rapidly rising: this period of time corresponds to the relaxation 
from the initial state (where H = 0) towards the metastable state. In the latter, ng(t) is almost constant for at least 
one, or even several decades of time. Then a decay of these plateau values sets in, which is due to the fact that too 
many much larger clusters have grown, the volume fraction of the system that is still in the metastable phase shrinks, 
and so less clusters of intermediate size (as studied in Fig. HOj) are observed. This behavior is qualitatively similar to 
previous studies (e.g. 0, [33[) which were based on the geometrical cluster definition, however. We also remark that 
in the case of Fig. ITUb the lifetime of the plateau extends up to t « 1000, one decade only, as expected from Fig. [j^i, 
since this case corresponds to a "lifetime" of the metastable state of only tms ~ 3000, and it is clear that only times 
t distinctly less than tms should be analyzed. 

While some of the previous work on the studies of the kinetics of cluster growth in metastable Ising models (e.g. 
[9l. l33l]) tried to use directly ni{t) to extract information on the validity of nucleation theory concepts, we here try to 
implement a different concept. Namely, we follow the trajectories of individual (large) clusters with respect to their 
size in time, {ti(t)} — > {l\{t + At)} where i is an index to label individual clusters. To ensure that the i'th 

cluster at time t + At is actually a descendent of the i'th cluster at the time f, we have to choose At small enough, 
and also record the location (center of gravity Xi{t)) and the components of the gyration radius 



RiAt) = < JJjj—J [ E <* - W) (^ «M) a I > - (27) 

where Xk, a is the a'th Cartesian coordinate for the fc'th lattice site belonging to the cluster with label i, consisting 
of £i(t) lattice sites. In each time step in which an analysis of the clusters is performed, the set of coordinates 
{Xi(t), Ri(t)} is recorded. 

Note that there occurs the difficulty that the number of large clusters is not constant during the simulation: clusters 
form and decay or split into parts, and since we know that f geom exceeds £ for each cluster, and the assignment of the 
"active bonds" according to Eq. ([5]) to identify from the geometrical cluster the associate physical clusters is a random 
process, some random shift of Xi(t) would occur even if we carry out two successive cluster identifications from the 
same spin configuration (At = 0). Of course, such shifts should be small in comparison with Ri(t), but as At is chosen 
nonzero it is clear that useful results are only obtained if At is small enough, and £ is large in comparison to clusters 
that correspond to typical thermal fluctuations [H, |6l| . Hence only clusters for which I > £ m i n are considered (for 
the temperature ksT/J = 3.0 we chose arbitrarily £ m j n = 10). So if by such criteria (for details see 89]) it is ensured 
that the z'th cluster with size £\{t + At) is a descendent of the i'th cluster with size £j(t) at time t, we can define a 
reaction rate T(£) as 

m _ ( stt^jm. ^ (28 ) 



1/2 
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Here the index £i(t) stands for an average over a sampling of all cluster trajectories recorded in the simulation, and 
a smoothing procedure of the (otherwise too noisy) data with a triangular smoothing function [89| was applied. Of 
course, in order to collect statistically significant data on T(£), it is necessary to perform many runs for each parameter 
combination (T, H, Hi) that is studied. We observe that the lifetime of the metastable stale in such runs is fluctuating 
dramatically, and so it is necessary to choose the run time of each run individually, rather than the same for all runs. 
It was decided to stop each run automatically when the largest cluster size ^™ ax (i) = L 3 /20 was reached. Of course, 
then only clusters with £ <C L 3 /20 could be studied, to avoid artifacts caused by this cutoff. 

(a) (b) 
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FIG. 11: (Color online) Cluster reaction rate T(£) versus cluster size £ for a bulk system at ksT/J = 3.0 without walls (a) and 
a system with walls at ksT/J = 3.0 and surface field Hi/ J — 0.1 (b). In each case several choices of the bulk field H are shown, 
increasing from bottom to top, for clusters grow more likely with larger field strength. Full curves are based on Eq. (|28[) , while 
broken curves are based on the approximation based on the use of the largest cluster only. In the heterogeneous case (b), both 
methods agree very well at lower fields because there is only one larger cluster in the system at a time. 



While Eq. (f2"5f is based on using all clusters £i(t) > £ m ; n at each time t, one can simplify matters by restricting 
the analysis only to the trajectory of the biggest cluster in the system [8^|. When one does this, one ignores possible 
problems from the fact that from time to time the identity of the largest cluster changes. Fig. [ITJ shows now typical 
results for T(£), using both this latter approximation and the method based on Eq. {28]) . Both methods yield similar 
trends, although they differ somewhat in detail (particularly in the case of homogeneous nuclcation in the bulk). 

What one expects theoretically for T(£) is a monotonous increase of T(£) with £, where T(£) is negative for clusters 
smaller than the critical cluster size £* while T(£) is positive for £ > £*. The method yields a second positive part for 
£ slightly larger than £ m i n . This is an artifact of ignoring clusters smaller than £ m ; n , which vanishes for £ m [ n — >• 1 [89| . 

If we identify the critical cluster size with the zero crossing of T(£) at the right hand side, and convert to the cluster 
volume V* according to Eq. (jTTj) . we obtain the data shown in Fig. [T5] It is seen that for a given value of the field 
H (or — /i cocx , respectively) classical nucleation theory underestimates the volume of the critical cluster: in other 
words, a given volume V* leads to a larger value of H . Since (according to classical nucleation theory and Eq. fl|)) 
the barriers scale as H~ 2 , this means when one studies nucleation barriers as function of cluster volume V* or cluster 
radius R*, one finds lower barriers in the simulation rather than predicted. Qualitatively, the data from the present 
analysis of cluster kinetics confirm the findings from the lever rule method of Winter et al. [55] , as far as homogeneous 
nucleation is concerned (Fig. \TTn ) , although some questions on systematic errors in both methods have not been fully 
settled. Nevertheless, the qualitative agreement between these quite different approaches is satisfactory. For the case 
of heterogeneous nucleation, however, for a given supersaturation the critical cluster volume predicted from cluster 
kinetics (Fig. ITT]) is distinctly larger than the corresponding results from the static methods. We have no explanation 
for this discrepancy. 
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FIG. 12: (Color online) Log-log plot of the critical cluster volume V* against the normalized chemical potential difference A/i 
for ksT / J = 3 for homogeneous nucleation in the bulk (a) and heterogeneous nucleation at a wall with several choices of the 
surface field Hi/ J, as indicated. The broken straight line is the prediction of the classical nucleation theory, Eq. ([3}, amended 
by the enhancement factor (1.064 at fcsT/J = 3) taken from Fig. [8] The dotted lines are the corresponding results based on 
the leverrule. Results for the method based on Eq. (|28[1 are shown in both plots, the method based on the biggest cluster is 
only drawn in (a), as both methods yield the same results in the heterogeneous case. In case (a) the lever rule data were taken 
from a system at linear dimension L = 15, and thus possibly affected by some finite size effects. Note also that in (b) the 
prediction of the classical nucleation theory and the leverrule results lie on top of each other. 
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IV. CONCLUDING DISCUSSION 

In the present work, we have studied aspects of nucleation theory by simulation of clusters and their dynamics, using 
the Ising (lattice gas) model on the simple cubic lattice. Both homogeneous nucleation and heterogeneous nucleation 
at planar walls (where a "surface field" acts) have been considered. 

Although many aspects of this problem have been studied before in works of various groups extending over several 
decades, most of the previous work is inconclusive since it relied on the use of the "geometric" cluster definition. 
We have given evidence that this geometric cluster definition does not yield correct results for large clusters at the 
temperatures far above the roughening transition temperature where the clusters have spherical shape; at temper- 
atures below the roughening transition temperature the geometric clusters and the "physical clusters" are basically 
indistinguishable, but due to the pronounced anisotropy effects a simple analysis of nucleation phenomena is not 
possible. 

However, in the limit of large droplet volumes V — > oo, where one can neglect any corrections to the decomposition 
of the droplet formation free energy into the bulk term plus a surface correction, one can compute the nucleation free 
energy barrier AF* from measuring the excess chemical potential A/i that is in equilibrium with a given V. Fig. [5] 
shows the enhancement of AF* with respect to the standard result for spherical droplets (using Eqs. ([3]), (0J). We 
thus show that in the Ising (lattice gas) model this enhancement gradually rises from unity as the temperature is 
lowered from the critical temperature, reaches almost 10% at T/T c = 0.6, and rises steeply below the roughening 
temperature towards the low temperature limit 6/n 1.91 (Fig. [5]). This enhancement reflects the consequences of 
the anisotropy of the intcrfacial free energy, such as the gradual crossover of the droplet shape from a sphere to a 
cube (Fig. [HJ) - 

On the other hand, we demonstrate that physical clusters do give consistent results, at least in the bulk when 
one is concerned with homogeneous nucleation. We show that in the limit where the droplets get macroscopically 
large, they converge against the simple lever rule predictions. However, we do find an (unexpected) surface excess 
in the particle number of such clusters also in this case. We also demonstrate the validity of the relation between 
chemical potential (of the supersaturated vapor) and the droplet radius that classical nucleation theory predicts for 
large droplets near the critical temperature. We also give evidence that the droplet-vapor interface is broadened due 
to capillary waves; we remind the reader that mean-field type theories and density functional theories [IJ E3 cannot 
include such capillary wave effects (which also should give rise to a correction term on the droplet formation free 
energy, not yet included in Eq. (j2J). 

We would also like to stress that many of our considerations can be carried over to a study of clusters in d = 2 
dimensions, where a construction as in Fig. Q] also holds. However, we expect two distinctions: (i) The roughening 
transition temperature Tr is zero, so the crossover of droplet shape from the circle to the square occurs without any 
singularity even for arbitrarily large droplets, (ii) Percolation coincides with the critical point, but geometrical clusters 
still are too large, and to describe nucleation, physical clusters defined via Eq. ([3]) should also be used. Of course, 
it would be very desirable to carry these considerations over to nucleation in off-lattice models of fluids. However, a 
precise analogue of Eq. (J5)) is still not known, and hence other concepts to define physical clusters [58| need to be 
used, if one wishes to study nucleation near the critical point. 

In the second part we present a first study of the time evolution of the cluster population based on the "physical 
cluster" definition. However, due to the large computer resources needed for this study, only data at a single tem- 
perature (ksT/ J = 3.0) are presented. In order to allow a comparison of this part of the study with our results on 
static properties of critical droplets, as studied in the first part of the paper, we use a criterion to estimate the critical 
droplet size from the balance between droplet growth and shrinking processes. In the case of homogeneous nucleation, 
the results obtained in this way are roughly compatible with the results obtained from the static lever rule method. 
Studying droplet volumes in the range from 100 to 200, clear deviations from the classical nucleation theory are seen, 
which can be attributed to a decrease of the nucleation barrier due to fluctuation effects. However, in the case of 
heterogeneous nucleation, a rather large discrepancy between the results of the statics and dynamics of droplets is 
found. This discrepancy is not understood yet, and must be left as a challenging problem for the future. 
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